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Abstract 

Interdependent networks are ubiquitous in our society, ranging from infrastructure to economics, 
and the study of their cascading behaviors using percolation theory has attracted much attention 
in the recent years. To analyze the percolation phenomena of these systems, different mathematical 
frameworks have been proposed including generating functions, eigenvalues among some others. 
These different frameworks approach the phase transition behaviors from different angles, and have 
been very successful in shaping the different quantities of interest including critical threshold, size 
of the giant component, order of phase transition and the dynamics of cascading. These methods 
also vary in their mathematical complexity in dealing with interdependent networks that have 
additional complexity in terms of the correlation among different layers of networks or links. In 
this work, we review a particular approach of simple self-consistent probability equations, and 
illustrate that it can greatly simplify the mathematical analysis for systems ranging from single 
layer network to various different interdependent networks. We give an overview on the detailed 
framework to study the nature of the critical phase transition, value of the critical threshold and 
size of the giant component for these different systems. 
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I. INTRODUCTION 


Systems consisting of mnltiple inter-connected networks with different types of links have 
received enormons attention in the recent years m, dne to its nbiqnitons applications in 
complex systems. Snch networks appear in the literatnre known as interdependent networks 
or mnltiplex networks. Stndies have shown that interdependent networks show distinct per¬ 
colation/ phase transition behaviors from single networks. In particnlar, an interdependent 
network is more vnlnerable to random attacks [10]. As many real world infrastrnctnre net¬ 
works can be classihed into interdependent networks mHH], the nnderstanding of their 
robnstness carries great practical signihcance. 

In a network consisting of links and nodes, one of the most important qnantities nsed 
to analyze its robnstness is the size of the giant component, which is dehned as the largest 
set of nodes that are connected with each other. When a network is nnder attack, i.e. a 
fraction 1 — p of nodes (or links) are removed, the size of the largest clnster shrinks. Usnally 
its size is a hnite fraction of the total nnmber of nodes in the network, nnless more than 
a certain fraction 1 — of nodes are removed - then the largest clnster (as known as the 
giant component) disappears and all of the clnsters become negligibly small. This phase is 
associated with the disintegration of the network. Hence the size of the giant component 
serves as an order parameter that is very nsefnl in stndying the phase transition behaviors 
and the robustness of the network structure. 

One of the original works in [T] provided a precise and powerful analytical solution to 
the phase transition behaviors. In their mathematical analysis, recursive mapping was used 
to track the percolation process in each stage of cascading failures. In some systems where 
correlations exist in dependency links [3-[9l [151 |T6], this method could lead to very compli¬ 
cated formulations and not always easy to solve. To study different network constructions, 
some of the other studies used different methods to achieve a relatively simpler analytical 
framework. In particular, the works in [H El [T71I20] used self-consistent equations of the 
converging probabilities to have an alternative approach to analyze for the critical behaviors 
on certain types of interdependent networks. Some of these methods can be extended to 
other scenarios. 

In this paper, we illustrate the use of one particular technique based on self-consistent 
probabilities |S1 ITTHTU] . and demonstrate that it could be applied to a wide variety of different 
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interdependent networks with minimum simplicity through surveying the literature in this 
held. This method focuses on the recursive representation of two central quantities dehned 
as the probabilities of hnding a link/node in the giant component. It is able to give a set of 
straight forward self-consistent equations describing the percolation behaviors without going 
through the cascading process [1] , and also deal with many correlated systems with simpler 
mathematics formulations. 

First we will illustrate the framework through the example of the single layer network. 
Next we extend it to multi-layer networks without degree degree correlations. Following that 
we extend the analysis to more complicated scenarios of partially correlated networks and 
degree-degree correlated networks. More complications are added to the case when multiple 
dependency links per node is introduced together with correlations, as well as single network 
with different types of links, also known as multiplex networks. 


II. SINGLE LAYER NETWORK 

The classic site percolation problem in a random network gives rich phase transi¬ 

tion phenomena for various networks structures. In the simplest case, we consider a random 
network without any correlations, and its degree distribution P{k) fully captures its struc¬ 
tural property. We start by introducing a key quantity x in the system; This x will be 
similarly dehned throughout this work and plays a central role in the mathematical analy¬ 
sis. If we randomly choose a link from the network and travel along one direction of the link, 
there is a probability x it would reach the giant component of the network, and probability 
(1 — x) it will not. (See Figurefor illustration). 

Suppose we randomly choose a link, and hnd an arbitrary node u by following this link 
in an arbitrary direction. The probability that the node u has degree k' is 

P{k')k' P{k')k' 

For this node u to be part of the giant component, at least one of its other k — 1 out-going 
links (other than the link we hrst picked) leads to the giant component. By calculating this 
probability, we can write out the self-consistent equation for x: 


3 




( 2 ) 


X = 




where 1 — (1 — ^ is the probability that at least one of the other k 

lead to the giant component, and is the probability that the node 


— 1 links of node u 
u has degree k. 



FIG. 1: Dehnition of x. A link (red color) is chosen and a node u (in green) is fonnd. From 
the three ontgoing links (dashed lines) of node A, one of them leads to to the giant 
component (represented by the symbol cxd). Since at least one of the three ontgoing links 
of u leads to the giant component, the red link leads to the giant component. And we 
define the probability of finding snch a red link as x. 


Therefore, for a randomly chosen node u, the probability that it is in the giant component 
is eqnal to the probability that at least one of its k links leads to the giant component. Thns 
we have: 

= (3) 

k 

where 1 — (1 — x)^ is the probability that none of the k links of node u leads to the giant 
component, and P{k) is the probability that node u has degree k. It is worth noting that 
/i°° is also the normalized size of the giant component, i.e. the fraction of nodes in the giant 
component. The above eqnations exactly eqnal to the results obtained by M. E. J. Newman 

and et.al in [23] . 

In the network percolation problem, when a fraction of 1 — p nodes are randomly removed 
from the network [231121], i.e. there is a fraction of p node remaining, we could apply the 
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previous equations with slight modifications. Assuming that the links of the removed nodes 
are still present on the network, the probability that a randomly selected link leads to the 
giant component is the same as before, given by ' [1 “ (1 “ But since only 

a fraction p of the nodes remain in the network, by calculating the probability that the 
randomly chosen link does not lead to the giant component, the self-consistent equation of 
X in En ([^ becomes: 

where 1 — (1 — is the probability that at least one of the k — 1 outgoing links of node 
u leads to the giant component, and is the probability that u has degree fc, same as 
before. The additional variable p in front is due to the fact that only a fraction of p nodes 
remain in the network after removing 1 — p nodes. 

Similarly, the probability that a randomly selected node is in the giant component is: 


(5) 

k 

It is known that in a single network, we usually only have second order phase transitions, 
such that there is no giant component when p is smaller than a critical probability pc- 
Above the threshold giant component appears and its size increases continuously from 0 
with increasing p. This means when p —)■ p[^, we have x ^ 0 and /r —)■ 0. When x —?• 0, by 
taking Taylor expansion of En Q, we obtain: 


X = 




{k) 


which leads to 




{k) 


( 6 ) 


(7) 


{k{k-l)) 

For Erdos-Renyi network, En Q yields p^J = l/{k), in agreement with the known result. 
For scale free network with 7 < 3, diverges, thus we would obtain pc = 0, also in 
agreement with the known result. 

The above system is based on node percolation, in which nodes are randomly removed 
until the giant component disintegrates. An alternative scenario is link (bond) percolation, 
in which links are randomly removed from the network. In this case, we still have the same 
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definition for x, and its equation remains the same as En because a randomly selected 
link has probability p to still remain in the network after removing a fraction of 1 — p links. 
The only difference is in En. in which we need to remove p on the right side: 

= ( 8 ) 

k 

This is due to the fact that all of the nodes remain in the network in link percolation, 
unlike the cases of node percolation that only a fraction of p remains. Hence we would 
obtain the same p[^ value for both node and link percolation, but different /i°° values. 

Before we proceed to interdependent networks, it is worth mentioning that other than 
second order phase transition mentioned above, there could also be hrst order phase transi¬ 
tion, and the critical threshold can be labelled as p^. For such phase transitions,when p < pi 
the size of the largest cluster is 0, and abruptly jumps to a non-zero value at p = p^. We 
shall see more of such examples later. 


III. MULTI-LAYER INTERDEPENDENT NETWORK 

A. Two Layer Interdependent Network 

In the original work of Ref pQ , generating functions was used to study the phase transitions 
in the two layer interdependent network. The system consists of two networks A and B, 
with degree distributions Pa(^) and Psik) respectively. Both networks A and B have N 
nodes, and each node in A is linked with exactly one node in B by a dependency link, 
and vice versa. The dependency link is different from the connectivity links within each 
network, in the way that once a node on one end of the dependency link is removed, the 
other node on the other network is also removed. This corresponds to the case where the 
failure of a power plant in the grid network will render the connected computer system to 
shut down due to the unavailability of electricity. Also, any node outside the giant cluster 
of its own network would fail since it is disconnected with the majority of the other nodes. 
In the dehned mutually connected giant component (MCGC), every node is in the giant 
component connected via the connectivity links in its own network, and its dependent node 
is in the giant component of the other network as well. Thus the MCGC is a steady state 
of the remaining network, such that no further cascading of failures would happen. 
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Here we present a simple method to study the phase transition behaviors using the 
formulation extended from the previous section. Following the dehnition in En ([^, we 
dehne x as the probability that a randomly chosen link in network A leads to the giant 
component. Analogously the probability that a randomly chosen link in network B leads to 
the giant component is y. 



FIG. 2: Dehnition of x and y in interdependent network. Two networks (A and B) both 
with N nodes are connected by dependency links, with one-to-one matching. Here node u 
is connected with u' via an interdependent link. A link (red solid line) in the network A is 
chosen and a node u (in green) is found following the link. Out of the two outgoing 
connectivity links (black dashed lines) of u, one leads to the giant component. Node u' in 
the lower network is connected with u via the dependency link (solid brown line), and it is 
also connected with the giant component via a connectivity link in network B. Since both 
node u and u' are connected to the giant component, we can be sure that the initial red link 
leads to the mutually connected giant component (MCGC). In network A, we dehne the 
probability of hnding a connectivity link leading to the MGGG as x. Similarly in network 
B, we dehne the probability of hnding a connectivity link leading to the MGGG as y. 

If a randomly chosen link in A leads to a node with degree A;, the node is in the MGGG 
only if at least one of its other k — 1 links lead to the giant component, and its dependent 
node in network B is also in the MGGG. Otherwise, this link will be not be in the mutually 
connected giant component, and be eventually deleted according to Ref [1] (see Figure 
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for a detailed illustration). Therefore by calculating the probability that a randomly chosen 
link in A leads to the MCGC, we would obtain: 


* = E - (1 - ■ E - (1 - (9) 

where “ (1 “ xY~^] is the probability that at least one of node u’s other k — 

1 connectivity links in network A leads to MCGC, and PB{k')[l — (1 — y)’^'] is the 
probability that at least one of the k' connectivity links of the dependent node u' in network 
B leads to the MCGC. 

Similarly we would have the probability that a randomly chosen link in B leads to the 
MCGC: 

V = E - (1 - - (1 - (“) 

k ^ k' 

Consequently, the probability that a randomly chosen node (either in network A or B) is 
in the MCGC is: 

= ■ E - (1 - (11) 

k k' 

which again is the normalized size of the mutually giant component. Note that we do not 
distinguish this value in different networks, because there is a one-to-one matching between 
nodes in A and B, so that /i°° is identical for both networks. 

When we randomly remove 1 —p fraction of nodes from network A, there is only p fraction 
of nodes left in A. Hence out of the original probability x that a randomly selected link leads 
to the MCGC, only a fraction of p nodes are actually remaining. It is easy to write down 
the new expression for x as: 

*=p ■ E ur^li - (1 - ■ E ^®(*=')ii - (1 - p)‘t (12) 

k ^ k' 

Analogously, the equation for y is 

p=p ■ E ^wii - (1 - p)‘"‘i ■ E ^■<(*^')ii - (1 - ^)‘t (13) 

t (*"») K 

At last, we arrive at the equation of /i°°, which is the probability that a randomly selected 
node in A (or B) is in the MCGC: 
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( 14 ) 


= P ■ - (1 - ^>‘1 ■ E - (1 - »)‘'l- 

k k' 

which is also the normalized size of the MCGC. 


In principle, Eqs. h^and 13 can be transformed into 


X = Fi(p,y), 


(15) 


and 

y = F2{p,x). (16) 

If we cannot get the explicit formula as above, the numerical computation always can be 
employed succesfully. 

Usually, the phase transition for the above system is of first order at the critical point 
(example given below). Therefore, at p = p^, the two functions x = Fi(p^,p) and 
y = F 2 {pi,x) meet tangentially with each other: 


dFi{pi,y) dF2{pi,x) 


= 1 . 


(17) 


dy dx 

For the first order phase transition, at the the critical point p^, the giant component is not 


0, implying that we cannot employ Taylor expansion to simplify Eqs. 1^ and 1^, but have 
to solve the polynomial equations directly. It could be very difficult to obtain the explicit 
formula for p^ except the most simple distributions, but numerical methods are possible. 


Example with random regular network 


For a simple example, we assume both network A and B are random regular networks 
with Ta(3) = Pb{^) = 1- It means every node in both networks have degree 3, and the 


nodes are randomly connected. The above Eqs. 12, 13 and 14 then becomes: 


X = p[l-{1- a;)^][l - (1 - p)^]. 

(18) 

p = p[l - (1 -p)^][l - (1 -x)^]. 

(19) 

“=p[l-(l-a;)=^][l-(l-p)3]. 

(20) 
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If a: 7 ^ 0 and 2 / 7 ^ 0, further simplification gives 


X = Fi{p,y) = 2- 

y = F 2 ip,x) = 2 


P[1 - (1-1/)^]’ 
1 

p[l — (1 — x)3] ’ 


( 21 ) 

( 22 ) 


Hence the requirement for of equation 17 can be written explicitly as 

3{l-yf 3(1-xf 


(23) 


Solving the above three equations gives us x = p ~ 0.5446, p^ ~ 0.7588 and consequently 
the mutual giant component size ~ 0.6329. 

The tangential requirement in equation 17 is presented in Figure When p = pi ^ 


0.7588, the curves from equations and 22 touch each other at x = p ~ 0.5446, where the 
slope of the two curves are equal; When p < pi, the two curves do not touch each other. 


and we only have the trivial solution of x = p = 0 from equations 1^ and IT This abrupt 
change in the size of the giant cluster corresponds to the hrst order transition, in which /r“ 
changes from 0 to 0.6329 abruptly at p = pi- This is illustrated through simulation results 
in Fig|^ 


B. N-layer interdependent network 

The case of N-layer interdependent networks o El EHl EMH] is an extension of the two 
layer scenario. Assuming there are N networks of equal number of nodes, and each node 
in a network is randomly connected with one and only one node in every other networks. 
Extending from En (|^, we obtain the probability that a randomly chosen link in network i 
leading the the MCGC as: 


*< = E _ (1 _ _ (1 _ ^ ,24) 

ki \ j^i kj 

where is the probability that a randomly chosen link in network i leading to node u 

has degree ki, [1 — (1 — Xj)*'*“^] is the probability that at least one of the other ki — 1 outgoing 
connectivity links of node u in network i leads to the MCGC, and Ylikj “ (1 “ 

is the probability that at least one of the kj connectivity links of node w in network j (node 
w and node u are connected by dependency link) leads to the MCGC. 


10 














FIG. 3: Solving for critical value of for random regular interdependent network with 
degree 3. The two curves represent Eqs. 21 and 22, At p = p], k. 0.7588, the curves from 


equations and [2^ touches each other at x = y ~ 0.5446, where the slope of the two 

curves are equal. 


Similarly extending En 11, we obtain the probability that a randomly selected node is in 
the MCGC as: 


/i~ = 


N 

( 25 ) 

n=l kn 

In the percolation problem, if 1 — p fraction of nodes are randomly removed from layer i, 


we can simply multiply En 24 and 25 by p, which is the fraction of nodes remaining in the 
layer after the attack: 


= p ■ E - (1 - ntE - (1 - V)‘'l} (26) 

ki ' j^i kj 

= Fi{p,xi,X 2 ,-■ ■,xn), (27) 

N 

=p - n{E^”(*=”)ii - (1 - ^”)‘"i}' (28) 

n=l k„ 

Example with Random Regular networks 
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p 


FIG. 4: First order phase transition on random regnlar interdependent network with 
degree 3. The red circles represent simulation result and the green line is the theoretical 


result from En 23 jumps from 0 to 0.6329 at p = = 0.7588, representing a first 

order phase transition. 


For a simple illustration, we let the networks to have the same degree distribution: 

F,(3) = P,(3) = 1, (29) 

i.e. every network is a random regular network with degree 3. 

Since the equations are symmetric, and there is a one-to-one matching between every 
node on each network, we have the relation 


Xi = X2 = ■ ■ ■ = Xn = X. 


(30) 


Now En 1^ simplifies into 


a; = p[l — (1 — a;)^][l — (1 — 


3iAf-l 


(31) 


Note that when iV = 2, we get exactly the En 18 with x = y. 


By bringing x to the right hand size, we could transform the En 31 into 


Fx = p[i - (1 - 2 :) ][i - (1 - xy] 


3lN-l 


X = 0. 


(32) 


12 












In this case, the critical value of pc can be understood as the smallest value of p such that 


En 32 has a real solution of x in the meaningful range of [0,1]. This means when p < p^, 


En 32 only has a trivial solution of a: = 0, when p > p^, there is more than one solution 
of a: > 0, and when p = Pc one unique solution of x exist . Thus at the critical point 
X = Xc,p = Pc, we would have the following relation fulhlled: 

dF^ 


dx 


r=Xc 0) 


(33) 


which leads to 


p{2{l - a:)[l - (1 - x)T" + 3(iV - 1)(1 - x)"[l - (1 - -1 = 0. (34) 

Solving the simultaneous equations and numerically, we are able to hnd out the value 
of pI- Notice that for any integer value of > 0, we have a solution of Pc and Xc in the 
range of [0,1], thus we always have a hrst order phase transition, but no second order one. 

IV. PERCOLATION ON MULTI-LAYER INTERDEPENDENT NETWORKS 
WITH DEGREE-DEGREE CORRELATIONS 


Usually in real-world networks, the connection through dependency links may not be 
random [3 El [151 ESI EHl jSSHSI]. In a general form, we can assume a joint probability 
P{kA, ks) for the dependency links to connect a node u with degree in network A and 
a node u' with degree in network B. In this case, we still assume that each node is 
connected with one and only one node in the other network through a dependency link. 

Instead of using the independent probabilities and Pb{P), fhe joint probability 


P{kA,kB) is used. Thus En 12, 13 and 14 become: 


X = p ■ 


E E - (1 - vf 


kA kg 


y = p'J2Yl ~~ [1 -(1 - 


kg kA 


X 


\^A] 




EE PikA, kB) [1 - (1 - xf^] [1 - (1 - y)'^^] 

kA kg 


(36) 


(36) 


(37) 
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In fact, En ^ ^ and are the more general representation of and ^ In the 

case of Ref [30], there is perfect correlation between the degrees of the two networks, i.e. 
P{kA,ks) = P{k) if kA = ks = k] else P{kA, ks) = 0. The above equations transform into: 


^ = P - E - (1 - P)‘l (38) 

p = p ■ E ^^(*^>11 - (1 - - (1 - (39) 


= P • ^(*=)|1 - (1 - 2')‘] [1 - (1 - P)‘l 


(40) 


A special case is when we have random regular networks for both A and B, and the 
results was discussed in the previous section since P{kA, ks) = P{kA) = Piks) = 1- 

The more general case of correlated systems of a multiplex network with different types 
of links were studied in Ref [8]. With similar argument, in the case of correlated N-layer 


interdependent networks, we could write down the equations of Xi and /i°° from En 27 and 


Xi=p- ^ P{ki,k2,...)^[l-{l-Xif^ 

kiM,- ' j¥-i 

= Fi{p,Xi,X2 ,-■-.Xn), (42) 

N 

p“=p' E (^^3) 

Ref [8] provided a general mathematical tool to solve for the critical points by using the 
Jacobian of the equations: 

det[J - I] = 0, (44) 


where J is the Jacobian matrix with = dFi/dxj. Solving En 42 and 44 gives the critical 
point value of and Xc for each layer of network in the system. 


V. PERCOLATION ON TWO LAYER PARTIALLY INTERDEPENDENT NET¬ 
WORKS 

In certain interdependent networks, not every node has a dependency link. It is more 
realistic to assume only a fraction of nodes from each network to have dependency links 
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j. And in such systems, both hrst and second order phase transitions may occur 
depending on the details of the networks’ structural properties. 

Let us assume two networks A and B with degree distribution PA^k) and Psik), and 
only a fraction of q nodes from each network is connected to nodes in the other network 
with dependency networks. For simplicity, we let each node to be connected with at most 
one other node through an dependency link. 

In order for a randomly selected link in A to lead to the MCGC, it must satisfy two 
conditions. First the node u it directly attaches to must have at least one of its outgoing 
connectivity links leading to the MCGC, and the probability is the same as the case of full 
dependency links given by “ (1 “ xY~^]. Secondly, for the case of network B, 

there are two scenarios: there is a probability 1 — q that node u is not connected with any 
node in B, then u is in the MCGC; there is a probability q that u is connected with a node 
u' in B, then at least one of the connectivity links of u' must also lead to the MCGC, and 


the probability for this is q PB{k')[l — (1 — y)’^']. Therefore the original En 12 becomes: 


x=p-Y^ 

k 


PA{k)k 

(k) 


[1 - (1 - xf-^] ■ {(1 - g) + g Pfyfc')[1 - (1 - yf]}. (45) 


k' 


The parameter p on the right hand side takes into account that after removal of 1 — p nodes 
in network A in the beginning of the attack, only a fraction of p node remain. 

It is worth noting that the calculation of y is not symmetric with x, because there is 
no one-to-one matching between a node in A and a node in B. For a node u' in B, the 
difference is in the case when it has a dependency node u m. A (probability g), at least one 
of the connectivity links of u must lead to the MCGC (probability PA{k')\l — (1 — x)^']). 


and u must not have been removed (probability p). Thus En 13 becomes: 


y= E wiwii - (1 - ■ {(1 - 9) + mY.PA(m - (i - («) 

There is no additional parameter p in front of the right hand side since we are not removing 
nodes from network B in the beginning. 

Again, due to the lack of symmetry in this case, the sizes of the MCGC in A and B are 
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expressed differently: 

A? = P • E - (1 - • {(1 - 9) + 9E - (1 - yt']} 

k k' 

p? = E Poimi - (1 - p)‘] ■ {(1 - 9 )+p9E - (1 - ^>‘'1} 


(47) 

(48) 


k' 


1. Example with random regular networks 


For a simple illustration, we use random regular networks for both A and B, and Pa{^) = 
-Pb( 3) = 1. Equations 45 and [4^ simplifies into 


X = p[2x — — g + g[l — (1 — 2/)^]}, 

V = pa - p^]{l -q + pqll - (1 -l)’]}. 


(49) 

(50) 


Note that, the cascading dynamics is not symmetric for network A and B for we only attack 
network A. For example, when g = 0, there is no dependency link between A and B] If we 
remove all of the nodes in network A, none of the nodes in B is affected and giant component 
still exist in B. Without the loss of generality, we study the case where 1 — p nodes are 
removed from A, and focus on the phase transition behavior in A. 

For the second order phase transition of network A, x = 0 at the critical point p = p^J. 
This means when p —)■ we have x —)■ 0. Note that This does not imply p —)• 0 (we use po 


denote this non-zero solution at the critical point). When x —)■ 0, from En 50 we have 

1 

2/0 = 2 - 


1 - g 


(51) 


Only the largest solution of po in fh® range [0,1] is the realistic solution, thus g must be in 
the range [0,0.5] for second order phase transition to occur. Usually for the more general 
cases, we could numerically get the solution of po by iterative calculations starting from a 
value close to 1. Note however, this solution only depends on g. Submitting po to Eq. 
and ignoring x^ when x —)■ 0, we have 


1 = Pc^2{l - g + g[l - (1 - go)^]}- 

Thus we can obtain the critical point value of p for second order phase transition: 

P" ^ 


211-9(1-(*)!])■ 


(52) 


(53) 
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For the first order phase transition of network A, again we cannot assnme x —)■ 0 as in 


the case of second order phase transition. Instead, by transforming Eqs. 45 and 46 we have 

1 


and 


X = Fi(p,y) = 2 


y = F 2 {p,x) = 2 


p{l-q + q[l - (1 -|/)3]}’ 


I — q + pq[l — (1 — 

As argned earlier, the critical point value satishes tangential requirement: 

d Flippy) dF2ipi,x) 


dy 


dx 


which can be written explicitly as 

3q{y - 1)2 


= 1 , 


3piq(x - If 


(54) 


(55) 


(56) 


pl{q[{y -1)3 +1] - g + 1 }^ {piq[{x - i)^ + 1 ] - g + 1 }^ 


= 1 


(57) 


Solving Eqs. 49, 50 and 57 for x, y and pf we are able to get the threshold value 
numerically. 

For any given value of g, which determines the fraction of dependency between networks 
A and B, we would have either a hrst order phase transition with critical threshold pf or 
a second order phase transition with critical threshold pf. Figure shows the plot of p^ 
and pf v.s. the change in g. The two curves of p^ and pf intersect at g = 0.4. Hence the 
percolation behavior is separated into two regions: When 0 < g < 0.4, it is second order 
transition; When 0.4 < g < 1, it is hrst order transition. 


2. Example with ER networks 

For a more general example, we suppose the degree distributions of network A and B 
are both poisson with average degree (fc/i) and {ks). In this case we can use the generating 
function formulation |23] to simplify the expressions. 

For networks A and H, the corresponding generating functions are Gf{x) = ^ PA{k)x^ = 
g(fcA>(a:-i) G^(p) = 'ffPB{k)y^ = The generating functions relating to 

the branching process are Gf{x) = gng[ G'f(p) = = 

g(*:s>(3/-i)_ Note that the generating functions of degree distribution and the the branching 
process are the same for poison degree distributions. 
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FIG. 5: and the giant component size Pa (.Pc) of network A at the critical point. 

We have and Pa{Pc) = 0 when q = 0.4. Pa{Pc) > 0 and increases with q when 

q > 0.4. It means that q = 0.4 is the bonndary of the second and hrst order phase 
transition. There is second order transition when q < 0.4 and hrst order transition when 
q > 0.4. q = 0.4 is the bonndary point for the system. 


Tims, the eqnations 45 and 46 can be written as 


X 


= p(l - e"'*-*’*)!! -? + ?(!- e"****')], 


(58) 


and 


y = (1 -g + pg(l -e-**-**")]. 


(59) 


Again, since to the natnre of poisson degree distribntions, the generating fnnctions of degree 
distribntion and the the branching process are the same, we have = x and = y. 
However, for other degree distribntions this relation is generally invalid. 


En 58 and 59 lead to the the explicit formnla 


(T = Fi(p,p) =-7^1og(l + ^ ^ 


{^a) 


pq 


pg(l — 


), 


(60) 
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and 


y = F 2 {p,x) = 


1 , .1 _ X 

{ks) Q P?(l — 


(61) 


For the second order phase transition of network A, a; —>■ 0 at the critical point. Using 

l/o = (l-g)(l-e-<'=«>"°). (62) 


Eq. 59 we have 


Snbmitting j/o to Eq. we can obtain the explicitly formnla of the second order phase 
transition critical point 

(63) 


P" = 




Here we nse the hrst order term in the Taylor expansion of 

For the hrst order phase transition of network A, nsing the tangential attachment of Eqs. 
[60] and (Ml we have 


dFi{pl,y) dF2{pi,x 


dy 


dx 


= 1 . 


(64) 


Again, Eqs. 58, 59 and 64 allow us obtain p^ numerically. Similar to the previous example, 


we can End out the value of q such that 


Pc = Pc ■ 


(65) 


This would allow us to End the boundary between Erst and second order phase transitions, 
i.e. the triple point value. 


Usually, the above three equations Ens. 60, 60 and 64 have no explicit formula allow 
us to solve it directly. In order to detect p{, we could brut-search for p G [0,1] according 


to following calculations. For a given p, we run Ens. and iteratively and obtain the 
solutions (Exed point). Then put this Exed point x, y into En. 64 If En. equal to 1, this 
p should be p[. For many more general case, such as the degree distribution is scale free, we 
have no explicit formula for En. then you have to use numerical way to degree partial 
derivative at the Exed point x, y. 

Numerical simulation could help us to End the critical points without solving the equa¬ 
tions [33]. As demonstrated in Egure|^ for 2nd order phase transitions, the second largest 
cluster size p 2 is maximum at the critical point; repeated simulation for diEerent value of p 
can be carried out to End out the peak /i 2 to identify pc. For 1st order phase transitions, 
the number of iterations (NOI) is at maximum; Thus one can identify the Pc as the point 


19 



















8 



P 


FIG. 6: Identification of critical point throngh simnlations. For 2nd order phase 
transitions (where q = 0.2), the maximnm size of the giant clnster appears at the critical 
point Pc- For 1st order phase transitions (where q = 0.8), the nnmber of cascading 
iterations is maximized at the critical point Pc- 

where maximnm NOI is located. Here iteration refers to the cascading of failnres from one 
network to the other. 

VI. SINGLE NETWORK WITH DIFFERENT TYPES OF LINKS 

Dependency links could also exist in single layer networks PESIEI], and the form of 
dependency could vary. This means, while the nodes in a single network are generally 
connected via connectivity links, some of the the nodes have mutual dependencies. 

Suppose for a given network, certain pairs of nodes are mutually dependent on each other 
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j. In this case, if node u is dependent on node w, w must and only depends on node u. 
Since the dependent nodes are in the same network, we could get the equation of x simply 


by replacing y with x in En 12 of the two interdependent networks.: 


X = 


( pE - (1 - - (1 - 


( 66 ) 


Note that the second p on the right side is there because both the node itself and its 
dependent node have the probability p to remain after the initial attack for they are in the 
same single network. Correspondingly, we can write down the size of the giant component 
as: 

= ( 67 ) 

For a more general case, a network could have dependency groups [31] - certain group 
of nodes have dependency relations (as shown in Fig. such that the removal of any one 
node would result in the removal of all the other nodes in the group. Given the probability 
distribution of the group size g{s) where s is the number of nodes in a group, we would 
obtain the following equation of x\ 


P = -tPE - (1 - P)*"']} ■ E - (1 - P)'"]}*-'- (68) 

Note that in this case, every node in the dependency group of size s needs to be in the 
giant component in order for the group to be in the giant cluster (else the whole dependency 
group would be removed). Here is the probability that a node is a group of size s, and 
— (1 — is the probability that every other s — 1 nodes is in the giant 

cluster. Using similar arguments, we could write down the equation for /i°°: 


P” = EE^{pE^W|1-(1-p)‘I}‘- (69) 

From here, we could study the critical phase transition behaviors by solving the equations 
using the techniques in the previous sections. 


VII. PERCOLATION ON TWO LAYER WITH MANY CORRELATED LINKS 

In certain networks including brain networks, it is observed that the dependency links 
are not always one-to-one matching, but could be one-to-many, and extensive correlations 
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FIG. 7: Single network with dependency groups. In this network, nodes in the same color 
are in the same dependency group. For example, if the green node on the bottom left is 
removed from the network, both of the two green nodes on the right would be removed due 
to their dependency relationship with the bottom left green node. 

and dependencies between nodes exist. Ref [7] discovered that brain networks are wired in 
such a way that stability is maximized. 

To examine the critical phase transition behaviors of such networks, additional parameters 
need to be dehned for specifying the structure. Assuming two networks A and R, we 
denote (k^) as the degree of connectivity links of a node in A{B), and k^^^ as its 

dependency degree that lead to nodes in the other network. The joint probability that a 
node M in A has k^ connectivity links and dependency links is denoted by PA{k^, ^^t)- 
The conditional probability PAB{k^\k^) is that given a node m in A with degree k^, the 
probability that any of its dependent node tc in R is of degree k^. Similar dehnitions carry 
over to network B for PB{k^, and PAB{k^\k!^). During an attack, a fraction of 1 — pi 
nodes in A and 1 — P 2 nodes in R are removed. This is in contrary to our previous examples 
in which only one network is being attacked at the beginning. 

Here we let xa(xb) be the probability that on following an arbitrary connectivity link 
in network A{B), we reach a node leading to the giant component. For an dependency 
link between node u with connectivity degree k^ and node w in B, ip^A is dehned as the 
probability that this link from u leads to the giant component, and yi.B is dehned similarly 
for a dependency link from R to A. Thus we have the following self-consistent equations: 
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( 70 ) 


Pi- 


E{ 

hA uA 
^in^'^out 


Xb=P2- X] { 

uB uB 


PAjkj, k^t)kt 

(kfn) 

PBjkfn, k!^^t)kg 

(kg) 


[1 - (1 - ^] ■ [1 - (1 - VkA 


k±- 


[1 - (1 - xbP" ^] ■ [1 - (1 - y^B 


(71) 


„ /^A j^A \^A 

where ^ "'JLaT "" is the probability that a randomly chosen link in network A leading to 

' in! 

node u has connectivity degree kg and dependency degree kg^^, [1 — (1 — XAg^"~g is the 
probability that at least one of the other kg — 1 ontgoing connectivity links of node u leads 
to the giant component, and [1 — (1 — yj^A is the probability that at least one of the 
kg^ dependency links of u leads to the giant component. 

For the probability yj^A and , we have 


y^t = P^(^PAB{kg\kgg)[l - (1 - XBf^}, ( 72 ) 

hB 

^in 

- (1 - XA^^}, ( 73 ) 

JuA 

^in 

where [1 — (1 — xa^^”-] is the probability that at least one of the kg^ dependency neighbors 
(in network B) of node u (which has degree kg) is in the giant component. Finally we can 
write down the probabilities that a randomly chosen node is in the giant component: 


Pa = Pi{ PA{kg^, kg^g[l - (1 - xa )^-] • [1 - (1 - ( 74 ) 

uA uA 

Pb = P2{ Y ^Bikfn, kgut)[^ - (1 - ■ [1 - (1 - ( 75 ) 

k-B uB 
^in^'^out 

which are straight forward. 

In general, since the above system is extremely complicated to have analytical solutions, 
numerical methods are usually preferred. Here we illustrate a simple yet efficient numerical 
method called binary search to detect the critical point for network A. The same method 
can be easily applied to network B. First we setup the initial starting points with pc- = 0 
and Pc+ = 1 , and let pc = 2 ^°^ . If Pa{Pc) = 0, we change the value of Pc- by letting 

Pc- = Pc] otherwise, we change pc+ by letting pc+ = Pc- A pc value with high precision could 
usually be reached with 20 such iterations as this algorithm converges exponentially. 
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VIII. CONCLUSION 


In this work, we have provided a specihc mathematical framework to review the criti¬ 
cal phase transition behavior of interdependent networks, otherwise known as Network of 
Networks (NON). Starting from single random networks, we have shown that by dehning 
two key mathematical qnantities - the probabilities of hnding a link/node in the final giant 
component, - one is able to directly write down the sets of self-consistent equations of these 
quantities without going through the iterative process of cascading failures in stages. This 
methodology greatly simplihes the mathematical analysis in complicated network structures, 
especially in very complex systems involving correlations and multiple dependency links per 
node. 

There has been many other works we have not included here. For example. Ref [9] has 
analyzed multiplex directed networks in the context of social networks. [351137] have studies 
evolutionary games. Interdependent networks with spatial constraint [38l [39] have been 
shown to exhibit unique phase transition behaviors, though we have not discussed them due 
to their analytical difficulties. In this work, our focus is to provide an mathematical overview 
using this specihc technique of simplihed self-consistent probabilities. The recursive mapping 
method [I] yields the same results, but we demonstrated that this particular method could 
greatly simplify the mathematical derivations for a wide range of complicated systems. 

Although this method proves to be applicable to a wide range of networks systems in 
studying their percolation behaviors, caution must be taken when implementing it. It is 
crucial that the self-consistent equations need to be carefully constructed, such that every 
component of the equations strictly follow the branching process underlying the percolation 
behaviors. 
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